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Summary. The kinetic theory of gases, including Granular Gases, is based on the 
Boltzmann equation. Many properties of the gas, from the characteristics of the ve- 
locity distribution function to the transport coefflcients may be expressed in terms of 
functions of the collision integral which we call kinetic integrals. Although the eval- 
uation of these functions is conceptually straightforward, technically it is frequently 
rather cumbersome. We report here a method for the analytical evaluation of kinetic 
integrals using computer algebra. We apply this method for the computation of some 
properties of Granular Gases, ranging from the moments of the velocity distribution 
function to the transport coefficients. For their technical complexity most of these 
quantities cannot be computed manually. 

1 Introduction 

The collision of elastic hard spheres of initial velocities vi and V2 is described 
by the collision law 



where v^^^ denote their velocities after the collision, vi2 = vi — V2 and 
e = fi2 1 r\2 with r\2 = ^ ^2 at the instant of the collision. For dissipa- 
tively colliding particles the coefficient of restitution e € (0, 1) characterizes 
the loss of relative normal velocity of the particles; e — \ represents elastic 
collisions, £ = stands for sticky collisions where the relative velocity af- 
ter a collision drops to zero. Consequently, the collision rule for dissipatively 
colliding particles reads 



v[ = ui - (ui2 • e ) e , 



V2 = V2 + (wi2 • e ) e 



(1) 




l+£ 

2 



{vx2 ■ e)e, 



V2 =V2 + 



l + £ 
2 



{vi2 ■ e)e. 



(2) 



The kinetic theory uses the concepts of direct and inverse collisions: In a direct 
collision, as described above, the particle velocities change by {ui,U2} — > 
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{w/i ^^2}- inverse collision is a collision which results in the velocities Vi/2, 
i.e. {vl',V2} -> {vi,V2}. 

The starting point of the kinetic theory of gases of elastic particles [1, 2] as 
well as dissipatively colliding particles (Granular Gases) [3] is the Boltzmann 
equation. It describes the evolution of the distribution function / (r, v, t) of 
the number of gas particles at time t in the infinitesimal phase-space volume 
dr, dv located at r, v. For dissipative gases it reads (e.g. [3-5]) 

^ + i/i • / (i7i, f) = a^'-i j dv2 J dee (-«i2 • el \vi2 ■ ^ X 

X [Xf{vi':t)f{v^',t)-f{v,,t)f{v2,t)] 

(3) 

which defines the collision integral I{f,f). Here a is the particle diameter 
and d is the dimension of the system. The integration is performed over the 
velocity V2 of the particle 2 of a colliding pair {1 — 2} and over the unit 
vector e. 0(x) denotes the Heaviside step-function. The Boltzmann equation 
describes the change of the population in the phase space due to the fiux of 
particles against the density gradient (second term in the Ihs of Eq. (3)) and 
due to collisions (rhs of Eq. (3)). 

The factor x in the gain term of the collision integral is specific for Granular 
Gases. It appears as a result of the transformation from the velocities u", Wj" of 
the inverse collision to the velocities vi,V2- Frequently it is assumed e = const, 
(we will call this simplified collision model, see section 3), the factor x is also 
a constant for this model, x = The case of molecular gases is obtained 

in the elastic limit e = 1. 

The collision integral has an important property: Given 

(V'(<)> = / dv,i^{v)f{v,t) (4) 

is the average value of some function ij: (v) and Aip (vi) = [ip (v- ) — -0 (t/^)] 
denotes the change of V' (vi) due to a direct collision. Then [1, 2] 

^(V'(i))= / dv,^{v,)^J{vi,t)^ I dvi^{vi)l{fj) 



dvidv2 J deO {-V12 ■ e) \vi2 ■ e\ f {vi,t) f {v2,t) 

X A[ijiv,) + ib{v2)] . (5) 

Hence, the time derivative of any avarage quantity tp (v) is expressed in terms 
of a function of the collision integral. The integration is performed over two 
velocities vi, V2 and over the unit vector e. The integrand is the product 
of two distribution functions and of the expression A [ijj (vi) + ijj {v2)], which 
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depends on the above three veetors since the after-collision velocities are given 
by Eq. (2). 

We will show below that all important properties of dilute homogeneous 
and non-homogeneous Granular Gases may be formulated in terms of expres- 
sions of the collision integral of the same structure as in Eq. (5). We will call 
these functions kinetic integrals. The structures of the kinetic integrals do not 
depend on the collision model, which allows to elaborate a general approach 
for gases of elastic particles, for dissipative gases with a simplified coUisional 
model e = const., as well as for gases of particles which collide via a realistic 
collision model as discussed in section 3. 

The aim of this paper is to describe a general approach for the analytical 
evaluation of kinetic integrals by means of symbolic algebra. To illustrate its 
application we compute various characteristics of homogeneous and inhomo- 
geneous Granular Gases. 



2 A simple example 

We consider a homogeneous Granular Gas where the velocity distribution 
function is independent of r. The granular temperature T is defined by the 
second moment of the velocity distribution function: 

^nTit)= Jdv^fiv^t), (6) 

which evolves due to [4-7] 

^ = -CT with C^V"S/^M2, (7) 

at a y m 

where C is the cooling coefficient of the homogeneous gas, n is the number 
density, m is the particle mass and the coefficient ^2 is the second moment of 
the dimensionless collision integral, with the general definition [4, 5, 7] 

/.jp = -i J dci j dc2 J deO {-C12 ■ e)\ci2 ■ e\ f (ci) f {C2) A{c^ + c^) . 

, (8) 

Here / (c) is the dimensionless distribution function 

fM^^jf^) , (9) 

which depends on the reduced velocity c = v/vt, with 



vr{t)^\l^ (10) 



being the thermal velocity of the gas. 
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For Granular Gases the velocity distribution function deviates from the 
Maxwell distribution. This deviation is characterized by the coefficients of the 
Sonine polynomial expansion: 



/(c) = 0(0) 



1 + ^ apSp (c^) 



(11) 



where the leading term (f> (c) = tt"''/^ exp (— c^) is the reduced Maxwell dis- 
tribution, ttp are numerical coefficients and the Sonine polynomials read for 
general dimension d: 



50 (x) = 1 

51 [x) = - 

etc. 



i(d + 2)x + i(d + 2)d 



(12) 



For small dissipation (e ^ 1) the deviations of the velocity distribution func- 
tion from the Maxwell distribution are small. Therefore, we approximate 



/(c)~(/.(c) [l + a2S2{^)\ 



(13) 



where we take into account that ai = due to the definition of temperature 
[4,5,7]. 

If we further neglect terms O (a|) , the second moment of the dimcnsionless 
collision integral reads 



dci / dc2 / deO {-Ci2 ■ e) \ci2 ■ e \ (j) (ci) (02) x 



X {1 + a2 [^2 (c?) + 52 {cl)] ]A{c\+ c 



This expression is of the form 

^2 / dci I dc2 / deO (-C12 • e ) \ci2 ■ e\ (j) (ci) (02) Expr(. . . ) 



(14) 



(15) 



with Expr = -i {1 + 02 [^2 {cl) + S2 (4)] }A{cl+ cl) 

We apply the following Maple program to obtain /i2 and, consequently, 
the cooling rate of the Granular Gas due to Eq. (7): 
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> restart ; 

> libname : =libname , ' /home/kies/KineticIntegral/ ' : 

> with(KineticIntegral) ; 

[DefDimension, DefJ, DefS, Klinit, getjexpr, unDefJ] 

> Klinit 0; 

> DefDimension (3) ; 

Dimension = 3 

Sonine polynomials defined 

> DefJO; 

Basic Integrals will be evaluated 

> Expr:=Delta2*(-l/2*(l+a2*(S(2,clp)+S(2,c2p)))) : 

> mu2 : =get Jexpr (0,Expr ,0) ; 

expressing 16 terms in form of Basic Integrals 

■- —L^^{£- 1) (1 + e) (16 + 3 a2) 



This example is intended to demonstrate the apphcation of symbohc algebra 
using the package Kineticlntegral .m for the evaluation of kinetic integrals 
without explaining the details. The mathematical method as well as the ac- 
cording Maple package are subject of the next sections. 
Keeping terms O (a|) wc obtain the result 

= \/2^(l-£') (^1 + ^02) , (16) 

which yields with Eq. (7) the temperature evolution law (Half's law) 

To 

(l+t/To) 



T{t)^ ,^ .2 (17) 



with the characteristic time tq 
where the mean collision time is given by 



r,{t)-'^4V^a'nJ^. (19) 
V m 

This example shows that with few lines only we can derive the temperature 
decay of a Granular Gas. Since the first part of the program is identical in 
all applications, only the main part is given below in the text. The program 
above then reads 
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1 Def Dimension(3) ; 

2 Def JO; 

3 Expr:=Delta2*(-l/2**(l+a2*(S(2,clp)+S(2,c2p)))) : 

4 mu2 : =get Jexpr CO ,Expr , 0) ; 



3 Granular Gases of viscoelastic particles 

It has been shown experimentally, e.g. [8], as well as theoretically [9] that 
the coefficient of restitution of realistic particles is not a material constant 
but it depends on the impact velocity. For the case of viscoelastic interaction 
the coefficient of restitution can be derived from the detailed analysis of the 
interaction of colliding spheres. The force between viscoelastic spheres of radii 
Ri and Rj which deform by ^ = \ fi — fj\ — Ri — Rj at the rate ^ is given as a 
superposition of Hertz's contact force [10] and a viscous force [11]: 

F = ( + UVii] ^ + Ua^ , (20) 



3(1 - ;y2) 2 ^ ^ 2 

where Y , v arc respectively the Young modulus and the Poisson ratio and 
= RiRjl {Ri + Rj). The dissipative constant A is a function of the bulk 
viscosities of the particles material (see [11] for details). Integrating New- 
ton's equation of motion for a particle collision the coefficient of restitution 
is obtained as a function of the normal component of the impact velocity 
g=\e- ivi-V2)\ [9,12] 

e = 1 - CiAa^/^gi/.s + ^^2A2a4/5g2/5 ^ ^ ^ ^ (21) 
5 

The numerical coefficient Ci reads 

1 21/552/5/^(21/10)- ^ 

Using the definitions of the thermal velocity, Eq. (10), and the reduced 
velocity c = v/vt the impact-velocity dependent coefficient of restitution, Eq. 
(21), may be written in terms of the reduced velocity: 

e = 1 - C\S' (t) |ci2 • e 1^/' + ^CfS'^ (t) |ci2 • e\'^' , (23) 

where 

6'it)^Aa'/'mt)f''^s(^^y\ (24) 

with 5 = Aa'^/^T^^^° and with To being the initial temperature. We will need 
this form below. 
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The derivation of the kinetic properties of Granular Gases of viscoelastic 
particles is significantly more complicated than for the case of the simplified 
coefficient of restitution (e = const.), see e.g. [5, 13]. The main reason of these 
complications is the factor x Eq. (3) which is not a constant anymore 
but a complicated function of the impact velocity [5] . The elaboration of the 
kinetics of Granular Gases of viscoelastic particles becomes very technical 
and was the main motivation for the application of computational symbolic 
algebra as presented here. 

In this paper we consider systems of three-dimensional spheres, i.e., for 
the case of viscoelastic spheres the particles interact by the force (20). Hence, 
for d = 2 we describe a gas of three-dimensional viscoelastic spheres whose 
motion is restricted to a plane. The coefficient of restitution for d = 2, i.e. for 
viscoelastic disks, has been derived recently [14, 15]. 



4 Evaluation of kinetic integrals 
4.1 Definition of the basic integrals 

The evaluation of kinetic integrals requires the computation of a sum of sim- 
pler integrals which we call basic integrals. These basic integrals may be com- 
puted analytically. Kinetic integrals are then decomposed into basic integrals. 
Although this decomposition is straightforward, the number of basic integrals 
may be as large as several thousands. Thus in practice it is not feasible to 
evaluate kinetic integrals manually. Let us consider a simple kinetic integral 
and show how it can be expressed in terms of basic integrals. In what fol- 
lows we elaborate the computation of kinetic integrals for general dimension 
d. Below the application of the method will be illustrated only for the most 
important case d = 3. The corresponding results for d = 2 can be obtained 
by changing the command Def Dimension (3) into Def Dimension (2) in the 
beginning of the programs. 

We assume that the dissipation is small enough to keep only first nonvan- 
ishing term in the Sonine polynomials expansion Eq. (11), i.e. we assume that 
the approximation 

fic)^q^ic)[l + a2S2{c^)]. (25) 

is adequate. With this approximation the moments of the dimensionlcss col- 
lision integral /ip reads 

fip^-^ J dci J dc2 j dee{-ci2 -6)1012 ■e\(f>{ci)(j) {c2)x 

X { 1 + 02 [52 (cl) + S2 (cl)] + al S2 (c^) 52 {cl) } A [cl + c|) . (26) 

To evaluate the kinetic integral in Eq. (26) it is convenient to use the center 
of mass velocity C and the relative velocity C12: 
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Cl^C+ ici2 , C2 = C - ici2 . (27) 

It is easy to check that the Jacobian of the transformation Eq. (27) equals 
unity. The product of the Maxwell distributions cj) {ci) (j) (C2) transforms into a 
product of a Maxwell distributions for the center of mass (an effective particle 
of double mass) and of a Maxwell distribution for the relative motion (an 
effective particle of half mass) 



1 / 1 o \ / 2 



d/2 

(j){ci)(t>{c2) ^ ; expl-^c^J (-) cxp(-2C2) 



= <^(ci2)<^(Cj . 
The terms in squared brackets in Eq. (26) hence transform into 

[52 (C?) + ^2 {cl)] = 

C'+{C- CI2) ' + ^4 + ic^cL - - Icj, + ^ . (29) 

To find the quantities Z\ (c^ + c^) {p = 2,4) we rewrite the collision law Eq. 
(2) in terms of the C and C12 too: 



c[ =C + ici2 - ^ (1 + e) (ci2 • e ) e 
C2 =C - ici2 + ^ (1 + e) (ci2 • e ) e . 



(30) 



Then straightforward algebra yields 



A{cl + cl) = -l{l-e')ic,,-ef (31) 



and 

Z\ {cj + 4) =2 (1 + ef (ci2 • e)' (c • e)' + ^ (l - e^)^ (c^^ . e)" 

e') (cl2 • e )^ c,l -{l-s^)C' (ci2 • e )^ (32) 
- 4 (1 + e) (C • C12) (C • e ) (ci2 • e ) . 

Combining the terms in Eqs. (29), (31), and (32) the integrands for fi2 and 
/X4 in Eq. (26) contain products of the Maxwell distributions for C and C12 
and the factors C*^, €^2, (c • C12)™ with exponents k, I, to, and similar other 
factors. Hence we define the basic integral by 

Jk,i,7n,n,p,a = Jdcu J dC J deO {~ci2 ■ e) \ci2 ■ e\^'^" 

X (C12) (C) ^^-42 {C ■ C12) " (C • e ) " (ci2 • ef , (33) 
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and express via these quantities. In particular /i2 reads in linear approxi- 
mation with respect to a2' 



f^2 = l{l-e' 



^0,0,0,0,2,0 + 0,2 [ «/4, 0,0, 0,2,0 + ^0,0,2,0,2,0 + «^0,4,0,0,2,0 



1 5 15 

+ - ^2,2,0,0,2,0 ^ 5J2, 0,0, 0,2,0 — ^ '/o, 2, 0,0, 2,0 + ^ >/o, 0,0, 0,2,0 



(34) 



The corresponding expression for ^4 contains already 35 basic integrals. Keep- 
ing the second order terms with respect to 02 yields expressions which contain 
60 basic integrals for ^2 and 328 basic integrals for fi4. To handle such ex- 
pressions one should obviously use computer algebra. 



4.2 Computation of the basic integrals 
Solution of the angular integrals 

The solution of the basic integral as defined by Eq. (33) is straightforward 
for n = 0, however, for n = 1,2 it requires some tricks (see, e.g., [16]). The 
presented approach may be generalised to n > 2, here we need only n = 0, 1, 2. 

In the d-dimensional space the Cartesian coordinates of a vector r = 
(xi, . . . , Xd) are related with the spherical coordinates, expressed by the ab- 
solute value of the vector r > and d — 1 angles 9i, 62, ... Od-i by the trans- 
formation 

Xi = r sin 01 sin02 ■ ■ • sin 0^-1 

d-l 

= rcos^rii-i JJ^ sin6'fc , m = 2,3,...d— 1 (35) 

k—m 

Xd = rcos0d-i , 

where 

< 6*1 < 27r 

- - (36) 
0<6'™<7r, m = 2,3, ...d-1. 

Hence, the infinitesimal solid angle reads 

d-l 

de^Y[sui''-Hkdek. (37) 

k=l 



Correspondingly we obtain for the surface area of a d-dimensional unit sphere 
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de^ dOi sin6l2(i6'2 / sir? e^de^ . . . j sm'^^^ ed-idOd-i 
Jo Jo Jo Jo 

r(|)r(i) r(2)r(i) r{^)r{i) 
rm r(2) ra) ■■ r'^^ 



2tt 



-'5 



d-2 , 



r(i) r(f)' 

(38) 



where we take into account that F {\) = 1, F(l/2) = ^Jtt and we use the 
trigonometric integrals 

7r/2 -, p n_ + l\ p tn+l\ 

sin' e cos" 6'd6' = -^r^ ■ (39) 







2 r(^ + i) 



and 



sin' = ^ 2 ;^ _ (40) 
^(5 + 1) 

Before calculating the Basic Integral we evaluate some more simple ex- 
pressions. We will need the angular integral 

/3,„ = j dee (e-i)(e-if = {-ir j deO (-e • (e • I") , (41) 

where vectors with a hat denote unit vectors (e.g. g = g/g). In Eq. (41) we 
change the variables e ^ — e to get the expression for j3m which will be mainly 
used in the calculations. Let us choose the coordinate axis OXd along g, so 
that e-g = cosOd-i- Using Eqs. (40) (39) the above properties of the Gamma 
function and Eq. (37) for the infinitesimal solid angle de we obtain 

Pm= de0{e-g^ ' ^) 



27r TT TT 7t/2 

d9i J 811162(192 J si-a^e^dOs... j sin''^^ ^d-i cos™ 9d-id9d-i 

%(i)r(i) rll)r{^) ir(^)r(^) 



(42) 



r(f) r(2) ■■■2 r{^) 

r f m+l \ 

^ (rf-i)/2 -' I 2 ; 

where we take into account that the integration over 9d-i is to be performed 
in the interval [0, 7r/2] due to the factor • (fj. Note that f3m does not 
depend on the vector g. Similarly we obtain 
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J de (e • I)" = J de [o (e • l) + (-e • |)] (e • l)" = [1 + (-1)™] /3„ 

(43) 

where we take into account that e is the integration variable and that Pm 
is independent of g. This allows to change the notation for the integration 

variable and for g: e—^ g and g C. 

Reduction of the basic integrals 

Wc also need the integrals 

j dCC^c^iC)^ J\'^-^+'^ (^^^ ' e-^^^ dC J dd 
and 

d99''4>{g) 



poo / 1 \ '^/^ r r ( ^ttL\ 

= 1 '""^' [ij ^-''^'dg jdg^2^'^^^ = 2r,^, (45) 

where again we have used J dg = J dC ~ J de ^ Qd according to Eq. (38). 
For 71 = the definition of the basic integral, Eq. (33), reads (with g = C12) 

JkM,p,c. ^ JdgJ dC<p{g)<p{C)C''g\C ■ gY''K^,p{g) , (46) 

where we introduce the integral 

X„,p(5)^ /dA/|e•.gr(e•g)^ (47) 



with the short-hand notation dfl — deO (— e - g) \e ■ g\. Similarly, for n = 1 
one can write 

Jk,i,m,i,p^a = JdgJ dC<l>{g)cl>{C)C'g' (c • 5)'" (c • /„,p(.9)) , (48) 
with the vectorial integral 

4,(g)^ / d^le\e■gr{e■gr ■ (49) 
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Finally, for n — 2: 

JkM.,p,c. = JdgJ dC^ig)(PiC)C''g' (c • .g)™ C • Ho..p{9) ' C (50) 
where the dyad Ha, pig) is given by 

Ha,p{g) = jdjlee\e-g\"{e-gY (51) 
and where we take into account that according to the vector-dyadic product 

C-ee-C=(c-ej(e-&j^(c-ej^ . (52) 

With the quantities introduced above we can express the basic integral step 
by step. 

Solution of the auxiliary vectorial and tensorial integrals 

The integral Ka,p{g) reads then 

KcAg) = / dee (-e • g) |e • g 1"+^ (e • gf 

= (-i)p,p+"+i I dee r.-. r^ ^y-"-^' (53) 



with the definition Eq. (41) of /3,„. Consider now the vectorial integral 
Ia,p{g) = gG{g). The result of the integration over e must be a vector. Due to 
the symmetry of the problem this vector should be directed along the vector 
g since there is no other preferred direction. Thus we write 

IcAg) = 9G{g) . (54) 
The fimction G{g) may be found by multiplying both sides of Eq. (49) by g: 

9-Ic.,p{9)^9^G{g)^ j dSei^e- g)\e- gr+' (e- gr+' = Ka,p+iig) . 

(55) 

Therefore, G{g) = g~'^Ka.p+i{g) and thus with Eq. (53) for Ka.p{g) we obtain 

laAg) = i-ir+'f3p+a+2 g'+'-g. (56) 

The integrand in Eq. (51) is a dyad proportional to the dyad e e. Since inte- 
gration over e does not change the symmetry, we expect that the components 
of the resultant dyad Ha.pig) composed of the components of the vector 
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g, the only vectorial quantity which is contained. Since the result may also 
have a non- vectorial part, the following form of Ha,p{g) is suggested: 

H[g) = A[g)gg + B[g)g^U, (57) 

where J7 is a unit dyad (i.e. the unit diagonal matrix) and the functions A{g) 
and B{g) are to be determined. Multiplying H scalar ly from both sides by g 
we obtain the first equation for A{g) and B{g): 

g ■ H ■ g = A{g)g^ + B{g)g^ 

= J dee i-e- g) |e • g 1"+^ (e • gf+^ = K^,p+2{g) , (58) 

where we use 

g-gg-g^ia^Y^g^ and g-U-g^g-g^g^. (59) 

Then taking the trace of H (i.e. computing the sum of diagonal elements) we 
arrive at the second equation: 

Tr H^A{g)g^+dB{g)g^^ J deO {-e ■ g)\e ■ g^+^e ■ gf ^ K^^p{g) . 

(60) 

In Eq. (60) (with the summation convention) we take into account 

T^T^ gg = gtgi ^ g^ ^ TTU = Sii = d, Tr ee = eic; = = l . (61) 

Solving the set of equations (58), (60) for the functions A{g) and B{g) and 
using Eq. (53) for Ka,p{g) we obtain 



(_l)PgP+^-l 



(d/3p+a+3 — Pp+a+l) gg + iPp+a + l — P p+a+s) g'^U 



(62) 



Solution for the basic integrals 



With the results for the scalar integral Ka.pig), foi" the vectorial integral 
Ia,p{g) and for the tensorial integral Ha^p{g) we can evaluate the basic inte- 
grals. For the first basic integral with n = we find with Eq. (53) 

Jk,iM,c. = {-lYPp+o,+i I dg j dC0(g)(/>(C)C'^+'"5'+™+^+"+' {p ■ l)"' 



{-ifPp+^+, j dcc^+^ci>{c) j d5.9'+'"+^+"+V(5) j dg{^-c) n^' , 

(63) 
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where we write the vectors g and C in the form g — gg and C = CC and for 
the integration over dg use 



dgFi{g)F2 (^) = ^ g^-^F^[g)dg J d^F2 



9''-'F,ig)dg J = J '^^Fiig) J d^F^ [^) Q^' 



(64) 



(recall that / = f^a)- The integral over d^" in Eq. (63) is [1 + (-!)"']/?„, 
according to Eq. (43), while the integrals over dC and over dg are given by 
Eqs. (44) and (45), respectively. As the result we obtain the basic integral for 
n = 0: 



= {-If [1 + (-1)™] 2'+'"+f+"+i|2^-i/3p+„+i/3™7fc+„,7,+„+p+„+i . (65) 

The evaluation of the basic integral for n = I and n ~ 2 may be performed 
by substituting Eqs. (56) and (62) correspondingly into Eqs. (48) and (50). 
The integrals which are obtained by this substitute are precisely of the same 
structure as in Eq. (63). Hence, the calculations for n = 1 and n ~ 2 are 
analogous to the case n ~ 0. The final result reads for n ~ 1: 

Jk,i,m,i,p,c. = i-iK' [1 + (-1)"+^] 2'+'»+p+"+ir2-i 

(66) 

and for n = 2: 



Jk,i,m,2,p,c = i-ir [1 + 2'+™+P+"+i [(d- l)/2rf]-i 

X 7fc+m+2 Jl+m+p+a+1 [(rf/5p+Q+3 " Pp+a + l) Pm+2 

+ (/?p+Q + l - Pp+a+s) Pm] ■ (67) 

The constants /3m , 7m and J?^ which are contained in the basic integrals have 
been defined above. 

We gave the solution of the basic integrals for n < 2 with arbitrary other 
indexes. As illustrated below, these basic integrals are mainly used in the 
kinetic gas theory. The present approach, however, may be straightforwardly 
extended. For example, the integral for n = 3 may be written (similar to (51) 

in terms of the tensor Ha',p^\g) = J d^eiCjCk |e • gj" (e • gf . Then this tensor 
and finally the basic integral may be solved using the same symmetry-based 
approach. 
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5 Computational formula manipulation to evaluate 
kinetic integrals 

Kinetic integrals have tlic structure 



dci / dc2 / deO (-C12 • e ) \ci2 ■ e\ (p (ci) (j) (02) Expr (...) , (68) 



where Expr (...) specifies a mathematical expression in terms of the pre- 
and aftercoUision velocities, the Sonine polynomials and their coefficients and 
other variables (see below). The evaluation of such integrals by means of ba- 
sic integrals Jk.i,m,n,p,a is straightforward, although for practical applications 
even simple looking expressions may result in sums of thousands of terms. 
Since the algebra follows an algorithm, wc can use a formula manipulation 
program, such as Maple to evaluate kinetic integrals. 

Before giving a systematic description we wish to describe the example 
from section 2 in detail, where we computed the second moment of the collision 
integral /i2 in linear approximation with respect to 02. 

The kernel of the kinetic integral Eq. (14) reads 



The computation in the listing in section 2 is started with restart to 
clean the memory of the computer from earlier results, then with libname . . . 
we extend the library path of Maple and load the library (with. . .). The 
next line is an output produced by maple, indicating that the commands 
Def Dimension, Def J, Def S, Klinit, getJexpr, and unDef J can be used now. 
Here starts the evaluation of our integral: First we initialize the computation 
by Klinit () which defines some internal variables. Then the dimension of the 
problem is specified. Here we wish to evaluate the Integral in three dimensions. 
With the line Def J() we require that the Basic Integrals will be evaluated, 
i.e., in the final result we generate real numbers according to (65)- (67), e.g., 
0,2, 0,2,0 = 6V27r. The line Expr : = . . . defines the kernel of the integral due 
to Eq. (68). Finally mu2:=getJexpr(0,Expr,0) starts the computation. The 
arguments of get Jexpr () as well as the syntax of Expr are described in detail 
below. 

We explain the Maple package Kineticlntegral .m in detail. To this end 
we comment the program line by line, bold numbers refer to line numbers in 
the listing. The package as well as the examples can be obtained at [17]. The 
program uses internal variables whose names follow a regular scheme: clp, 
c2p, cla, and c2a stand for the pre-coUision and after-collision velocities of 
the particles. The unit vector e= E = r^/r 12 and the dimensionality DIM 
are further variables. Scalar products are named by xDOTy, where x and y are 
vectors. An example is clpD0Tc2p= ci • C2. A full list of the internal variables 
is given in the table below. All these symbols can be used to specify the kernel 
Expr of the kinetic integral. 





(69) 



i4o i . Foscnel, IN . 


V. Brilliantov 


clp 


C^ ' 
^1 • 


pre-collision velocity of particle 1 


c2p 




pre-coUision velocity of particle 2 


cla 


c[: 


after-collision velocity of particle 1 


c2a 




after-collision velocity of particle 2 


cl2 


Cl2 


relative velocity 


C 


C: 


centre of mass velocity 


cl2D0Te=eD0Tcl2 


Cl2 


• e 


CD0Tcl2=cl2D0TC 




Cl2 


clpDOTclp 




=(72+Cil/4 + C-Ci2 


c2pD0Tc2p 




- (72 + Ciy4 - C • Ci2 


clpD0Tc2p =c2pD0T 


clp C?L • 


C2 c.l/A 


claDDTcla 


(c1 


f =(72+Cil/4-(l-e2)/4(ci2-e 



+ C • ci2 - (1 + e) (C • e) (ci2 • e ) 

c2aDDTc2a {c!^f ^ C"^ + c^l /4 - (l - e^) /4 (ci2 • e )^ - 

- C • ci2 + (1 + e) (C • e) (ci2 • e ) 

claDDTc2a=c2aD0Tcla c[ ■ c!^ = & - c^^jA + (l " /4 (ci2 ■ 
clpDOTcla=claDOTclp ci • c{ = (7^ + Ci|/4 + C • Ci2- 

-(1 + e)/2 (C • e) (ci2 • e ) - (1 + £)/4 (ci2 • e 
c2pDDTc2a=c2aD0Tc2p C2 • = (7^ + Cj^|/4 - C • ci2 + 

+ (1 + e)/2 (C • e) (ci2 • e ) - (1 + £)/4 (ci2 • e f 
claDDTc2p=c2pD0Tcla c[-ci=&- 

-(1 + E)I2 (C ■ e) (ci2 • e ) + (1 + e)/4 (ci2 • e 
clpD0Tc2a=c2aD0Tclp cx ■ c'^ ^ & - c-^2/^+ 

+ (1 + e)/2 (C ■ e) (ci2 ■ e) + {1 + £)/4 (ci2 • e 

Delta2, Delta4 see Eqs. (31), (32) 

DIM dimension 

CI, dprime describes the coefficient of restitution of 

viscoelastic particles, see (23) !!!! 
abscl2D0Te5 \ci2 ■ e\^^^ 

The variables Delta2 and Delta4 are defined in the initialization KlinitO: 



Klinit :=procC) 

global Delta2,Delta4: 
unprotect('Delta2' , 'Delta4') : 
Delta2 : = (l/2*epsilon-2-l/2) *cl2D0Te-2 : 
Delta4 : = (-4*epsiloil-4) *CD0TG*cl2D0Te*CDDTcl2 
+ (l/8-l/4*Gpsiloil"2+l/8*Gpsiloil"4)*cl2DDTe~4 
+ (-l/4+l/4*epsilon-2)*cl2-2*cl2D0Te-2 
+{2*epsllon"2+2+4*epsllon)*cl2D0Te"2*CDDTe~2 
+(epsilon"2-l)*C"2*cl2D0Te"2: 
protect {'DGlta2' , 'Delta4') : 
end: 
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The kernel Expr may contain Sonine polynomials. They are defined by Def SO 
according to Eq. (12). 



12 


DefS:=proc() 




13 


global S: 




14 


if (assignedCDlM) ) then 




15 


unprotect ( ' S ' ) : 




16 


S:=proc(a: : integer, x: lalgebraic) 




17 


local r,xx: 




IS 


xx:=subs(clp=C+l/2*cl2,x) : 




19 


xx:=subs(c2p=C-l/2*cl2,xx) : 




20 


xx:=subs(cla=C+l/2*cl2-l/2*(l+epsilon)*cl2DDTe*e 


xx) : 


21 


xx:=subs(c2a=C-l/2*cl2+l/2*(l+epsilon) *cl2DDTe*e 


xx) : 


22 


xx:=expandCxx"2) : 




23 


xx:=algsubs(C*cl2=CD0Tcl2,xx) : 




24 


if (a=0)then 




25 


r:=l: 




26 


elif (a=l) then 




27 


r:=-xx+DIM/2: 




28 


elif (a=2) then 




29 


r:=xx-2/2 - (DlM+2)*xx/2 + (DlM+2) *DIM/8 ; 




30 


else 




31 


ERROR C ' invalid argument a) : 
fi: 




32 
33 


r:=factor(r) : 




34 


end: 




35 


protect ( ' S ' ) : 




36 


printf C "Sonine polynomials defined\n"); 




37 


else 




38 


ERROR( 'Dimension undefined') : 




39 


fi: 




40 


end: 





First at line 14 it is checked whether the dimensionality DIM is specified, oth- 
erwise Def S reports an error at line 38. The syntax of the Sonine polynomials 
is S (a, expr) , a specifies the order of the polynomial and expr is its algebraic 
argument which may contain the pre- and after-collision particle velocities 
clp, cla, c2p, and c2a. In lines 18-21 these velocities are translated into 
expressions of the center of mass velocity C, the relative velocity cl2, the 
coefficient of restitution epsilon, the unit vector e and the scalar product 
cl2D0Te, due to Eq. (30). There is an important difference in the notations 
of the Sonine coefficients due to Eq. (12) and in the Maple program: 

S (a, expr) = Sa (cxpr^) . (70) 

This difference of the notation is taken into account by line 22 where the entire 
term Expr with all the described substitutions is squared. In the squared ex- 
pressions there appears the term C •Ci2 which is again substituted by CD0Tcl2. 
The program evaluates Sonine polynomials of the order a = {0, 1,2} which 
is done in line 24-29. Any incorrect input a is acknowledged by an error 
message, line 31. 

All of the defined internal variables, including the Sonine polynomials are 
accessible within the Maple program. In general, it is of not much value to list 
them since according to the substitution of variables even small terms become 
impressive expressions. For debugging purposes, however, this feature may be 
useful. 
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There arc two difTercnt output modes: first we may be interested in the 
result of an integration in terms of basic integrals Jk,i,m,n,p,a a-nd second 
(usually) we are interested in obtaining real numbers. The translation of the 
basic integrals into real numbers is given by Eqs. (65)- (67). These equations 
are implemented in Def J(): 



Def J : =proc C) 
global J: 
unprotect C J) : 
if (assignedCDlM) ) then 
J : =proc (k: : integer , 1 : : integer ,m: : integer ,n: : integer , 
p : : integer , a : : integer) 
local JO, Jl, J2, beta, gamma, omega: 

beta:=m->Pi- ( (DIM-1) /2) *GAMMA{ (m+1) /2) /GAMMA ( (m+DIM) /2) : 
gamma : =m->2- (-m/2) *GAMMA ( (m+DIM) /2) /GAMMA (DIM/2) : 
omega : =2*Pi- (DlM/2) /GAMMA (DIM/2) : 

J0:=(k,l,m,p,a)->(-l)"p*(l+(-l)"m)*2~(l+m+p+a+l)/omega*beta(p+a+l) 

*beta(m) *gamma(k+m) *gamma(l+m+p+a+l) : 
Jl :=(k,l,m,p, a) ->(-l)"(p+l) *(!+(-!)" (m+1) )*2"(l+m+p+a+l)/omega*beta(p+a+2) 

*beta(m+l) *gamma(k+m+l) *gamma(l+m+p+a+l) : 
J2:=(k,l,m,p,a)->(-l)"p*(l+(-l)"m)*2"(l+m+p+a+l)/omega*gamma(k+m+2) 

*gamma(l+m+p+a+l)*( (DlM*beta(p+a+3)-beta(p+a+l) )*beta(m+2) 

+(beta(p+a+l)-beta(p+a+3))*beta(m))/ (DIM-1) : 
slmpllfy(l/2*J0(k,l,m,p,a/5)*(n-l)*(n-2)-Jl(k,l,m,p,a/5)*n*(n-2) 

+l/2*J2(k,l,m,p,a/5)*n*(n-l)) : 
end: 

protect ( ' J ' ) : 

printf ("Basic Integrals will be evaluated\n") : 
else 

ERROR( 'Dimension undefined'): 
fi: 
end: 



After calling def J all results are given as real numbers by evaluating the basic 
integrals, e.g., 



> DefJO; 

Basic Integrals will be evaluated 

> BIdemo:=J(0,0,l,l,l,l) ; 



11 2(3/5) Trta/^) 

Bldemo := 



100 1 N p/ 9 
sinl — -K) I I — ) 



> evalf (Bldemo) ; 



2.811423645 



The counterpart of DefJO is unDef J(). After calling unDef J any result is 
given in terms of basic integrals. 



unDef J : =proc () 
global J : 

if (assigned( J) ) then 

unprotect (J) : 

unassign( ' J ' ) ; 
fi; 

printf ( "Basic Integrals will not be evaluatedXn" ) : 
end: 
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In the described programs DefS and Def J we have already used the variable 
DIM which contains the dimensionality of the Basic Integral. This variable is 
set by Def DimensionO : 



75 


Def Dimension : =proc Cd : :integ 


er) 


76 


global DIM: 




77 


unprotect ( ' DIM ' ) : 




78 


DIM:=d: 




79 


printf ("Dimension = %d\n" 


,d): 


80 


protect('DIM') : 




81 


DefSO : 




82 


if CassignedC J) ) then 




83 


Def JO : 




84 


fl: 




85 


end: 





After changing the variable DIM in line 78 the Sonine polynomials which are 
specific for any dimension have to be redefined (line 81). For the case that 
Def J() has been called before, the basic integrals have to be redefined as well, 
hne 83. 

The heart of the Maple package is the program getJexprO. This pro- 
gram has three arguments: epstype specifies the type of the particle interac- 
tion. epstype=:0 stands for constant coefficient of restitution and epstype=l 
stands for viscoelastic particles. The second argument is the kernel of the 
integral expr which is a function of the particle velocities before and after 
the collision, the dimension, the Sonine polynomials and their arguments, the 
Sonine coefficients, and scalar products of the vectorial values. In fact expr 
may contain all variables which are listed in the table on page 146. The third 
argument describes the output mode of the program. It will be described later. 
The text of the program getJexprO is listed below: 



100 
101 
102 
103 
104 
105 
106 
107 
108 
109 
110 
111 
112 



get Jexpr : =proc Cepstype : : integer , expr : : algebraic , 
ou: : integer) 
global di,ddi: 

local Expr , termJ , it , term,ko ,epow, aa,k, 1 ,m,n,p , a2pow, 

vareps , Clexpl ,ClomegaO , C12omegal : 

Expr : =subs (c2pD0Tclp=clpD0Tc2p , expr) ; 

Expr :=subs(c2aD0Tcla=claD0Tc2a, Expr) ; 

Expr : =subs (clpDOTcla=claDOTclp , Expr) ; 

Expr : =subs (c2pD0Tc2a=c2aD0Tc2p , Expr) ; 

Expr : =subs (c2pD0Tcla=claD0Tc2p , Expr) ; 

Expr : =subs (c2aD0Tclp=clpD0Tc2a, Expr) ; 

Expr:=subs(clpD0Tclp=(C*C+l/4*cl2*cl2+CDDTcl2) ,Expr) ; 
Expr:=subs(c2pD0Tc2p=(C*C+l/4*cl2*cl2-CDDTcl2) ,Expr) ; 
Expr:=subs(clpD0Tc2p=(C*C-l/4*cl2*cl2) ,Expr) ; 

Expr : =subs (claDOTcla= (C*C+l/4*cl2*cl2- (l-epsilon*epsilon) /4*cl2D0Te*cl2D0Te+CD0Tcl2 

-(l+epsilon)*CDDTe*cl2DDTe) ,Expr) ; 
Expr : =subs (c2aD0Tc2a= (C*C+l/4*cl2*cl2- (l-epsilon*epsilon) /4*cl2D0Te*cl2D0Te-CD0Tcl2 

+(l+epsilon)*CDDTe*cl2DDTe) ,Expr) ; 
Expr:=subs(claD0Tc2a=(C*C-l/4*cl2*cl2+(l-epsilon*epsilon)/4*cl2DDTe*cl2D0Te) ,Expr) ; 
Expr : =subs (claDOTclp= (C*C+l/4*cl2*cl2+CDDTcl2- (l+epsilon)/2*CD0Te*cl2D0Te 

-(l+epsilon)/4*cl2DDTe*cl2DDTe) ,Expr) ; 
Expr : =subs (c2aD0Tc2p= (C*C+l/4*cl2*cl2-CDDTcl2+(l+epsilon) /2*CDDTe*cl2D0Te 

-(l+epsilon)/4*cl2DDTe*cl2DDTe) ,Expr) ; 
Expr : =subs (claD0Tc2p= (C*C-l/4*cl2*cl2- (1+epsilon) /2*CD0Te*cl2D0Te 

+(l+epsilon)/4*cl2D0Te*cl2DDTe) ,Expr) ; 
Expr : =subs (clpD0Tc2a= (C*C-l/4*cl2*cl2+( 1+epsilon) /2*CD0Te*cl2DDTe 

+(l+epsilon)/4*cl2D0Te*cl2DQTe) ,Expr) ; 



150 T. Poschel, N. V. Brilliantov 



130 


k 


=degre 


e(term,C) : 


131 


1 


=degre 


e(term,cl2) : 


132 


m 


=degre 


G(term,CD0Tcl2) 


133 


n 


=degre 


G(term,CDOTe) : 


134 


P 


=degree (term , cl2DDTe) 



if (epstype=0) then 

Expr : =expcmdCExpr) : 
elif (epstype=l) then 

vareps:=l-Cl*dprime*abscl2D0Te5+3/5*Cl~2*dprime~2*abscl2D0Te5'2: 

Expr : =expcindCalgsubs Cepsilon=vareps ,Expr) ) : 
else 

ERRDR( ' invalid argument' ,epstype) : 
fi: 

Expr :=expand(siniplifyCExpr)) : 
termJ : =0 : 

printf ("expressing %d terms in form of Basic Integrals\n" , 

nops (Expr) ) : 
for it from 1 to nops (Expr) do 

term : =opCit , Expr) : 

ko:=coeff sCterni,{epsilon,abscl2DQTe5,C,cl2,CD0Tcl2,CDQTe,cl2DDTe,a2>) : 
epow:=degreeCterni,epsilon) : 
aa:=degreeCterni,abscl2D0Te5) : 



a2pow:=degree(term,a2) : 
term J : =terniJ+epsilon~epow*ko*a2"a2pow* J (k , 1 , m , n ,p , aa) : 

di [it] :=ko*epsilon~epow*cl2DQTe"aa*C~k*cl2~l*CDDTcl2'm*CDDTe~n*cl2DaTe'p*a2"a2pow: 
ddi [it] :=terDi: 

if Csimplify(di[it]-ddi[it]) !=0) then 

ERRORC error evaluating term', it,', see di [] , ddi[]'): 
fi: 
od: 

termJ : =simplif y CtermJ) : 
if (ou=l) then 

Clexpl:=GAMMA(3/5)*sqrt(Pi)/C2-(l/5)*5-(2/5)*GAMMA(21/10)) ; 

termJ : =subs (Cl=Clexpl , termJ) : 

expcindCevalf CtermJ) ) : 
elif (ou=0) then 

factor Csimplify(termJ)) : 
elif (ou=2) then 

ClomegaO:=omegaO/(2+sqrtC2*Pi)*2-(l/10)*GAMMAC21/10)) : 

C12omegal : =omegal/ Csqrt (2*Pi) *2' C 1/5) *GAMMA (16/5) ) : 



taylor(termJ,dprime=0,3) : 
normal (simplify CtermJ , GAMMA, trig) ) ; 
subs (Cl~2=C12omegal , termJ) ; 
subs(Cl=ClomegaO,termJ) ; 
normal (simplify CtermJ) ) ; 



termJ 
termJ 
termJ 
termJ 
termJ 
else 

ERRDR( ' invalid argument ', ou) 
fi: 
end: 



First, in lines 91-112 the kernel of the integral expr is processed. All variables 
are expressed in terms of epsilon, C, cl2, CD0Tcl2, CDOTe, cl2D0Te, a2. If we 
deal with viscoelastic particles (epstype=l), the coefficient of restitution is a 
function of the relative velocity. It is expressed in our variables by Eq. (23), 
where the coefficients 5' describes the material properties of the viscoelastic 
particles and C\ is a numerical constant, see Eq. (21). Therefore, for this case 
e has to be replaced by an expression in our variables which is done in lines 
116-117. We introduce a new internal variable abscl2D0Te5= |ci2 ■ e |"^^^. 

The entire kernel of the integral is then in line 121 represented as a sum 
of polynomials in the variables C, cl2, CD0Tcl2, CDOTe, cl2D0Te, a2, epsilon 
(for e = const.), and abscl2D0Te5 (for viscoelastic particles). These are pre- 
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cisely the variables which define our basic integrals for which we know the 
solution (see Eqs. (65)-(67)). For any of these terms we determine in the loop 
(lines 125-142) the indices of the according basic integral. The terms are then 
successively picked out of the sum by line 126 and stored in term. The next 
line determines the coefficients of the variables in this term, they are stored 
temporarily in the list coeff s. In lines 128-134 the powers of the variables 
in the present term are determined, i.e., epow is the power of epsilon in 
I term I ,aa is the power of abscl2D0Te5, etc. These powers determine the in- 
dices of the basic integral. In this way all terms of the sum are transformed in 
line 136 into expressions of the basic integrals and the according prcfactors. 
By processing the terms in the loop the variable term J grows and finally it 
contains the full sum in terms of basic integrals. 

The vector di [] in line 137 is for testing purposes. It reconstructs the 
present term by using the powers of the variables as computed above and its 
prefactors. If the reconstructed term differs from the original term (line 139) 
an error has occurred. This may happen by an typing error for the kernel 
of the integral. The program then reports the error and the number of the 
concerned term. 

In line 143 the integral is completely expressed in terms of basic integrals. 
There are different output modes, specified by the third argument ou of the 
program getJexprO. If ou=0 the output is given in symbolic form, using the 
Gamma- function and its properties. The result is analytically exact, however, 
the output may become very large. For ou=l the result is given by a floating 
point number. Before, the numerical constant Ci is replaced in line 146 by 
its analytical expression Eq. (22). For the choice ou=2 the result is given 
as a Taylor expansion using the numerical constants a;i/2, which make the 
expression more compact. These constants read: 

ujo = 2V2^2^^^°r (^^^ Ci w 6.485 , ivi = V2^2^^^r (^y^ w 9.285 . 

(71) 

Finally the defined commands have to be bundled to a package and stored 
with the name Kineticlntegral .m. This is done by the last lines of the 
program: 



162 
163 
164 



168 
169 



Kineticlntegral : =table C) : 
Kineticlntegral [Klinit] : =eval (Klinit) : 
Kineticlntegral [Def J] : =eval (Def J) : 
Kineticlntegral [unDef J] : =eval (unDef J) : 
Kineticlntegral [getJexpr] : =eval (get Jexpr) : 
Kineticlntegral [Def Dimension] : =eval (Def Dimension) : 
Kineticlntegral [Def S] : =eval (Def S) : 
save Kineticlntegral, 'Kineticlntegral .m' ; 



The presented program generates the maple library Kineticlntegral. m. 
To evaluate integrals of the type Eq. (68) the library can be invoked and used 
as demonstrated in the example in section 2. 
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6 Kinetic integrals in the kinetic theory of Granular 
Gases 

6.1 Homogeneous cooling state 

Temperature decay and velocity distribution function for the 
simplified collision model (e = const.) 

In the first stage of the evolution of an initially homogeneous gas of dissipative 
particles the gas preserves its homogeneity, however, its temperature decays 
persistently. This state of the evolution is called the homogeneous cooling 
state. In this state the velocity distribution function does not depend on the 
space coordinate r. Hence, the homogeneous cooling state is characterized by 
the gas temperature T as defined by Eq. (6). 

In our introducing example in section 2 we have derived the evolution of 
temperature of a Granular Gas. The deviation of the velocity distribution 
function from the Maxwell distribution is characterized by the Sonine coeffi- 
cients Ui (see Eq. (11)). If the deviations from the Maxwell distributions are 
small only the first nonvanishing term in the expansion Eq. (11) contributes 
considerably to the distribution function, i.e. we use approximation (13). For 
the simplified collision model (e = const.) the second Sonine coefficient 02 is 
related with the moments of the collision integral by [4. 7] 

5Ai2(l + a2)-A'4 = 0. (72) 



This equation is solved by the program below: 



1 


Def DimensionC3) ; 




2 


Def JO ; 




3 


Expr:=Delta2*(-l/2*(l+a2*(S(2,clp)+S(2,c2p)) + a2-2*S (2 , clp) *S (2 


c2p))): 


4 


mu2:=getJexpr(0,Expr,0) ; 




5 


Expr:=Delta4*(-l/2*(l+a2*(S(2,clp)+S(2,c2p)) + a2-2*S (2 , clp) *S (2 


c2p))): 


6 


mu4:=getJexpr(0,Expr,0) ; 






A2 : =solve {5*mu2* ( l+a2) -mu4=0 , a2) : 




8 


AA2 : =unapply (Re (A2 [2] ) , epsilon) : 




9 


A2NE : =epsiloii-> (16* (l-epsilon) * (l-2*epsllon"2) ) 




10 


/(81-17*epsilon+30*epsilon"2*(l-epsilon) ) ; 




11 
12 
13 


plot ( [AA2 , A2NE] , . . 1 , linestyle= [1,4], 
legend= ["nonlin . theory" , "lin. theory (v. Noije, Ernst)"], 
font= [COURIER, 20] , thickness=2 , labels= [e , a2] , color=black) ; 





As a result we obtain the second Sonine coefficient as a function of the coef- 
ficient of restitution as drawn in Fig. 1. The full analytic expression for a2{e) 
can be found in [4] . 

Temperature decay and velocity distribution function for gases of 
viscoelastic particles 

For Granular Gases of viscoelastic particles under the adiabatic approximation 
[18] the cooling rate and the second Sonine coefficient are computed in the 
same way as for e = const. The according program reads 
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. 2 



0.15 



a2 



0.1- 



0.05- 




Fig. 1. Second Sonine coefficient 02 over the coefficien of restitution e. Full line: 
solution of Eq. (72) (nonlinear theory [4]), dashed line: linear solution of Eq. (72), 
i.e. only linear terms with respect to a2 are kept in this equation [7]. 



DefDimensionC3) : 
Def JO : 

Expr:=Delta2*(-l/2*(l+a2*{S(2,clp)+S(2,c2p)) + a2-2*S (2 , clp) *S (2 , c2p) ) ) : 
mu2 : =get Jexpr (1 ,Expr , 2) ; 

Expr:=Delta4*(-l/2*(l+a2*(S(2,clp)+S(2,c2p)) + a2-2*S (2 , clp) *S (2 , c2p) ) ) : 

mu4 : =get Jexpr (1 ,Expr , 2) : 

A2 : =solve (5*iiiu2* (l+a2) -mu4=0 , a2) ; 



which yields ^2 on line 5 in terms of 02 and then 02 on line 9. Note that on 
line 9 we need to take a physical root 02 oc (5'. The result for 02 reads 



0.2 = a,2\5' -f 022^' + 



021 



-0.388, 



022 



12063 27 



(73) 

2.752. (74) 



20V27r 640000 tt 40 ^2^ 

Correspondingly the cooling coefficient, expressed via //2 (see Eq. (7)) reads 



c 



(c.o'5'-t^25" + ...) , 



where 



W2 



500 



(75) 



(76) 



which yields for the temperature evolution in linear approximation with re- 
spect to (5 



m 

Tn 



t 

1 + — 

To 



-5/3 



(77) 



where we introduce the characteristic time 
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r,' = ^go<5r,(0)-i = ^q,64<j'nJ^, (78) 

O V Tfl 



with the mitial mean colhsion thiie Tc(0) and the constant 

<Zo = 2V5r(^|)^«. 0.173. (79) 

The evaluation of a2 requires already the solution of 786 basic integrals. The 
above results for the temperature evolution and second Sonine coefficient have 
been reported in [4, 5]. 

Coefficient of self-difFusion for the simplified collision model 
(e = const.) 

If we assume Molecular Chaos in the homogeneous cooling state the velocity 
correlation function 

{vmt')) - (vitr) cxp (-^^) (80) 

is related with the coefficient of self-diffusion 

D{t) = -T,{t) (81) 
m 

via the velocity relaxation time t„. This quantity reads [19] 

X J dcidc2 j de6> (-C12 ■ e) |ci2 ■ e I / (ci,t) /(c2,t) (1 + £) (ci2 ■ e )^ . 

(82) 

Again we recognize the kinetic integral in Eq. (82) which allows for the appli- 
cation of the Maple package. Using the program 



Def DimensionCS) ; 
unDef JO ; 

Expr : =tau_c_inv/ (12*sqrt (2*Pi) ) * (1+epsilon) *cl2D0Te*cl2DDTe* (l+a2*S (2 , clp) ) 

*(l+a2*+S(2,c2p))) : 
Def JO ; 

tau_v_inv:=getJexpr(0,Expr,0) ; 



we obtain the velocity relaxation time 

2 



.,.,)^<i±il(l + §) .-W (83) 
which yields with Eqs. (81) the self-diffusion coefficient 
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2D„{t) 



(84) 



(1 + e) (1 + ^02) 



2 ' 



where 




Tit)- 



(85) 



is the Enskog sclf-difFusion coefficient of the eiastic gas. 

The self-diffusion coefficient in the linear approximation with respect to 
02 fias been computed for the simplified collision model in [20]. 

CoefHcient of self-difFusion for gases of viscoelastic particles 

A similar analysis can be performed for gases of viscoelastic particles: 

1 DefDimensionC3) : 

2 DefJC): 

3 Expr: = (l+epsilon)*cl2D0Te*cl2D0Te*(H-a2*(S(2,clp)+S(2,c2p))+a2"2*S(2,clp)*S(2,c2p)) : 

4 hattauV:= get Jexpr CI ,Expr , 2) : 

5 tauVinv : =1/ (2*DIM) *sqrt (2*T/m) *n*sigma~ (DIM-1 ) *hatt auV : 

6 Expr:=Delta2*{-l/2*(l+a2*{S(2,clp)+S(2,c2p))+a2-2*S(2,clp)*S(2,c2p))) : 

7 mu2 : =get Jexpr (1 ,Expr , 2) : 

8 Expr:=Delta4*{-l/2*(l+a2*{S(2,clp)+S(2,c2p))+a2-2*S(2,clp)*S(2,c2p))) : 

9 mu4 : =get Jexpr (1 , Expr , 2) : 

10 tmpA2:={solve{mu4=lim2/DIM*4*DIM/4*(DIK+2)*(H-a2) ,a2)}: 

11 a2:=tmpA2[l] : 

12 a2 : =convert (a2 , polynom) : 

13 tauVinv : =convert Ctaylor (tauVinv ,dprime ,3) , polynom) : 

14 Dvis:=T/ni/tauVinv: 

15 DEnsk:=DIM*GAMMA (DIM/2) *sqrt(T/m) / (4*Pi- ( (DIM-1) /2) *n*sigma- (DIM-1) ) : 

16 simplify (convert (taylor (Dvis/DEnsk,dprime , 3) , polynom)) ; 



At line 14 we need to take a root 02 cx S'. The coefficient of self-diffusion for 
gases of viscoelastic particles in linear approximation with respect to 02 has 
been reportet in [19]. Here we give the complete solution with respect to 02 
which becomes feasible due to symbolic algebra: 



6.2 Inhomogeneous Granular Gases 
Transport coefficients 

So far we have considered Granular Gases in the homogeneous cooling state. 
To describe the non-homogeneous gas a number of perturbative approaches 
has been developed. The most successive arc the Chapman-Enskog approach 
[1, 2] and the Grad method [21]; here we consider the former one. 

In the Chapman-Enskog approach it is assumed that (after a short relax- 
ation time) the time and space dependence of the velocity distribution function 



Dit)=Do{t) 1+ W-c,o<5' + 




625537 Lo^ 
20480000 V 



1851 /Y 
14080 V^' 
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/(r, i/, t) occurs only via density n(r, t), flow velocity u(r, t) and temperature 
T(f,t). These quantities are, respectively, the zeroth, first and second order 
moments of the velocity distribution function; n(r, t) and u{r^ t) are expressed 
in terms of /(r, v, t) in the same way is temperature in (6) [1, 2]. These three 
fields evolve according to the hydrodynamic equations (e.g. [13, 18]), 

dn 

— + V • {nu) = 
at 

o — > 

+ W+(nm)"^V-P = (87) 

at 

where P denotes the pressure tensor and q is the vector of the heat flux. In the 
Navier-Stokes hydrodynamics the pressure tensor and the heat flux depend 
linearly on the field gradients: 



Pi] = p5ij - r] ViUj + Vjiii - -%V • u 



2 



with the hydrostatic pressure p = nT in the dilute limit. Here 77 is the vis- 
cosity coefficient, k is the coefficient of thermal conductivity and coefficient 
jjL characterises the heat flux due to the density gradients and does have an 
analogue for elastic gases. 

The transport coefficients ry, k, fi may be calculated using the Chapman- 
Enskog approach adopted for the dissipative gases. The most complete anal- 
ysis of this problem for the gas with a simplified coUisional model with a 
constant e have been performed by Brey and coauthors [22,23]. They have 
expressed transport coefficients for the dissipative gases in terms of functions 
of the collision integral. The generalization of this approach for the dissipative 
gases with the velocity-dependent restitution coefficient has been done in [18], 
where it has been shown that the transport coefficient may be expressed in 
terms of the same functions of the collision integral as for the case of e = const. 

Transport coefHcients for the simplified collision model 

Using the notations of [18] the viscosity coefficient reads for the case of e = 
const, (see also [22]), 



1 mT 1 



where ^2 has been found above (see Eq. (16)) and with the function of the 
collision integral 
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n^= J dci J dc2 J dee (-C12 • e ) |ci2 • e | /(°Hci)x 

X (/.(C2) A, (c2)Zi [A, (Cl) + A, (C2)] , (90) 

where Aj (c) is the traceless tensor, 

Dij (c ) EE c,Cj - i%c^ . (91) 

In Eq. (90) we recognize the kinetic integral. For the simphfied coUision 
model (e = const.) it is solved by the program 



Def DimensionC3) ; 
Def JO ; 

Expr :=(-l/DIM*c2pD0Tc2p*Delta2+claDDTc2p*claDDTc2p+c2aD0Tc2p*c2aD0Tc2p 

-c2pDDTc2p*c2pD0Tc2p-clpDDTc2p*clpDDTc2p) * (H-a2*S (2 , clp) ) : 
OmegaETA:=getJexpr(0,Expr,0) ; 



yielding 

Similarly the coefficient of thermal conductivity may be expressed as 

15 rT~ l + 2a2 



l2„ = -V27r(l + £)(3-e)(l-g) , (92) 



8^252(^7) V 2m ^2 + i^^K ' ^^^^ 
with the function 1?^ defined as 

= j dci J dc2 J de0(-ci2-e)|ci2-e|/(°)(ci)x 

X 0(c2)5 {C2)-A [S (ci) + S{c2)] , (94) 

and with the vectorial quantity, 

- f.2 5 



Sic) ^[c'--]c. (95) 



The coefficient /x may be also written in terms of i?^: 



The function f}^ is computed by the program 



Def DimensionC3) ; 
Def JO ; 

Expr : = (c2pDDTc2p- {DIM+2) /2) * (claD0Tc2p*claDDTcla+c2aD0Tc2p*c2aD0Tc2a 

-c2pD0Tc2p*c2pD0Tc2p-clpD0Tc2p*clpDQTclp) * (l+a2*S (2 , clp) ) : 
OmegaKAPPA:=getjGxpr(0,Expr,0) ; 



158 T. Poschel, N. V. Brilliantov 



We obtain 

/7^,. ^ /49 - 33£ 19 - 3e \ 
a = -\/2^(l + e) + ^^^2 j . (97) 

In linear approximation with respect to a2 the transport coefficients have been 
reported in [22]. Substituting Eqs. (92,97,16) into Eqs. (89,93,96) we obtain 
the gencraUzed result for the coefficients r;, k and /i, which accounts for the 
terms of the next order, O {a^. 

Transport coefficients for Granular Gases of viscoelastic particles 

The same function appears in the theory of Granular Gases of viscoelastic 
particles [18]. We evaluate this fmiction by 



Def DimensionC3) ; 
Def JO ; 

Expr:=(-l/DIM*c2pD0Tc2p*Delta2+claDDTc2p*claDDTc2p+c2aD0Tc2p*c2aDDTc2p 

-c2pD0Tc2p*c2pD0Tc2p-clpD0Tc2p*clpDQTc2p) * (l+a2*S (2 , clp) ) : 
OmegaETA : =get Jexpr ( 1 , Expr ,2); 



and obtain the corresponding coefficient 

= -{wo + 5'wi^5'^W2) (98) 



with 



Wo = 4\/27r ^1 - 

= "° ( 1^ " 550"^) 
/ 97 679 
V 165 44000 



Similarly, by 



DefDimensionC3) ; 
Def JO ; 

Expr : = (c2pD0Tc2p- (DIM+2) /2) * (claDDTc2p*claDDTcla+c2aD0Tc2p*c2aDDTc2a 

-c2pD0Tc2p*c2pD0Tc2p-clpD0Tc2p*clpDQTclp) * (l+a2*S (2 , clp) ) : 
0megaKAPPA:=getJexpr(l,Expr,2) ; 



we obtain the other function 

= - (uo + 5'ui - 5'^U2) (100) 

with 



uo = 4V2^ + ^02^ 



17 9 

ui = Wo ( — - 77;7;a2 ) (101) 



5 500 



U2 = OJl 



/1817 1113 \ 
Vlio" " 352000"v 
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For the gas of viscoclastic particles the viscosity coefficient rj is found as 
the solution of the equation [18] 



where C is given by (75). We solve this equation as a perturbative expansion 
with respect to the dissipative parameter S' , 



using the program 



1 Def DimensionC3) : 

2 Def JO: 

3 Expr:=Delta2*{-l/2*(l+a2*{S(2,clp)+S(2,c2p))+a2-2*S(2,clp)*S(2,c2p))) : 

4 mu2 : =get Jexpr (1 ,Expr , 2) : 

5 mu2 : =convert Cniu2 , polynom) : 

6 Expr:=Delta4*(-l/2*(l+a2*(S(2,clp)+S(2,c2p))+a2-2*S(2,clp)*S(2,c2p))) : 

7 mu4 : =get Jexpr (1 ,Expr , 2) : 

8 tmpA2:={solve{mu4=mu2/DIM*4*DIM/4*(DIM+2)*(H-a2) ,a2)}: 

9 a2:=tmpA2[l] : 

10 a2 : =coiivert (a2 , polynom) : 

11 xlO:=2*sigma"(DIM-l)/DIM*n*mu2*sqrt(2*T/m) : 

12 ExprETA:={-l/DIM*c2pDQTc2p*Delta2+claDQTc2p*claDQTc2p+c2aD0Tc2p*c2aDDTc2p 

13 -c2pD0Tc2p*c2pD0Tc2p-clpD0Tc2p*clpDQTc2p)*(l+a2*S(2,clp)) : 

14 DmegaETA : =get Jexpr ( 1 , ExprETA , 2) : 

15 DmegaETA :=coiivert(QmegaETA, polynom) : 

16 EqETA:=xi0*eta0*(l/2+3/5*dprime*etaH-7/10*dprime"2*eta2) 

17 +4/ (DIM-1) / (DIM+2) *etaO* (l+dprlme*etal 

IS +dprlme"2*eta2)*n*sigma"(DIM-l)*sqrt (2*T/m)*DmegaETA+n*Ti 

19 EqETA:=taylor(EqETA,dprlme,3) : 

20 etaO:=solve(coeff (EqETA,dprime,0) ,etaO) ; 

21 etal :=solve(coef f (EqETA , dprime , 1) ,etal) ; 

22 eta2:=solve(coeff (EqETA, dprime, 2) ,eta2) ; 

23 omegal :=9.28569: 

24 omega0:=6.48562: 

25 evalf (etal) ; 

26 evalf (eta2) ; 

27 etaVIS:=etaO*(l+dprime*etal+dprime~2*eta2) ; 

28 plot CetaVIS/ etaO , dprime=0 . .0.5, color=black , linestyle= [1] ) ; 



The result reads 



Similarly, the coefficient of thermal conductivity may be obtained from 
the equation [18], 




(102) 




(103) 




(104) 




(105) 



which is also solved perturbatively as an expansion 




(106) 



The result is found by the Maple program 
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Def DimensionC3) : 
Def JO : 

ExprKAPPA : = (c2pD0Tc2p- (DIM+2) IT) * (claDaTc2p*claDaTcla+c2aD0Tc2p*c2aDDTc2a 

-c2pD0Tc2p*c2pD0Tc2p-clpD0Tc2p*clpD0Tclp)*(l+a2*S(2,clp)) : 
OmegaKAPPA : =get Jexpr ( 1 , ExprKAPPA , 2) : 
OmegaKAPPA : =convert (OmegaKAPPA , polynom) : 

Expr:=Delta2*{-l/2*(l+a2*{S(2,clp)+S(2,c2p))+a2-2*S(2,clp)*S(2,c2p))) : 
mu2 : =get Jexpr (1 ,Expr , 2) : 

Expr:=Delta4*(-l/2*(l+a2*(S(2,clp)+S(2,c2p))+a2-2*S(2,clp)*S(2,c2p))) : 
mu4 : =get Jexpr (1 ,Expr , 2) : 
mu2 : =convert Cmu2 , polynom) : 

tmpA2 : ={solve (mu4=mu2/DIM*4*DIM/4* (DIM+2) * (l+a2) , a2) } : 

a2:=tmpA2[l] : 

a2 : =convert (a2 , polynom) : 

mu2 : =convert Cmu2 , polynom) : 

xlO : =2*sigma" (DIM-1) /DIM*n*mu2*sqrt (2*T/m) : 

xlOl : =coef f (xlO .dprlme , 1) / (2/DIM*sqrt (2*T/m) *sigma" (DIM-1) *n) : 
xl02 : =-coef f (xiO ,dprime ,2) / (2/DIM*sqrt (2*T/m) *sigma" (DIM-l) *n) : 

TdxlOdT : =2/DIM*sqrt (2*T/m) *sigma" (DIM-1) *n* (3/5*xi01*dprime-7/10*xi02*dprime"2) : 
a21 : =coef f (a2 , dprlme , 1) : 
a22:=coeff(a2, dprlme, 2) : 

Tda2dT:=l/10*a21*dprlme+l/5*a22*dprime"2: 

OmegaKAPPA :=convert(taylor (OmegaKAPPA, dprlme, 3) , polynom) : 

EqKAPPA : =DIM* (kappa*TdxlOdT+xlO*TdkappadT+xlO*kappa) 

+4/ (DIM+2) *kappa*n*slgma" (DIM-1) *sqrt (2*T/m) *OmegaKAPPA 

+ (DIM+2) *DIM/4*n*2*T/m* ( l+2*a2+Tda2dT) : 
kappa: =kappaO* (l+dprlme*kappal+dprime~2*kappa2) : 
TdkappadT:=kappa0*(l/2+3/5*dprime*kappal+7/10*dprime"2*kappa2) : 
EqKAPPA : =tayIor (EqKAPPA , dprlme , 3) : 
kappaO : =soIve (coeff (EqKAPPA , dprlme ,0) ,kappaO) ; 
kappal : =soIve (coeff (EqKAPPA , dprlme , 1) ,kappal) ; 
kappa2 : =soIve (coeff (EqKAPPA , dprlme , 2) ,kappa2) ; 
omegal :=9.28569: 
omega0:=6.48562: 
evalf (kappal) ; 
evalf (kappa2) ; 

plot(kappa/kappaO,dprime=0. .0. l,color=bIack,llnestyIe=[l] ) ; 



Thus we found the coefficient of thermal conductivity 
75 / T 

K= -\\ (l-|-0.393(5'4-4.904<5'2 + ...) . (107) 

64(7^ V 7rm ^ ' 

Completely analogously wc obtain the remaining transport coefficient [x: 




1.229(5' + 1.415(5'2 + ...) (108) 



7 Conclusion 

We have reported a method for the analytical evaluation of functions of the 
collision integral, which wc call kinetic integrals, as they are encountered in 
the kinetic theory of Granular Gases. This method is based on computational 
symbolic programming, which operates with a set of known basic integrals. To 
illustrate the application of this method we evaluate various properties of the 
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Granular Gases, ranging from the moments of the velocity distribution func- 
tion to the transport coefficients. Two different models of the Granular Gases 
are considered: The simplified model, where particles collide with a constant 
coefficient of restitution, and the more realistic model where the restitution 
coefficient depends on the impact velocity, as it follows from the collision 
law for viscoelastic particles. We demonstrate that evaluation of some kinetic 
integrals may not in practice be done manually, while the reported method 
allows to obtain analytical expressions for these quantities. In the examples 
we generalize some previously known results for the kinetic properties of the 
dissipative gases. 
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